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The  flooded  agglomerate  model  has  found  prolific  usage  in  modeling  the  oxygen  reduction  reaction  within 
the  cathode  catalyst  layer  of  a  polymer  electrolyte  membrane  fuel  cell  (PEMFC).  The  assumption  made 
in  this  model  is  that  the  ionomer-coated  carbon-platinum  agglomerate  is  spherical  in  shape  and  that 
the  spheres  are  non-overlapping.  This  assumption  is  convenient  because  the  governing  equations  lend 
themselves  to  closed-form  analytical  solution  when  a  spherical  shape  is  assumed.  In  reality,  micrographs 
of  the  catalyst  layer  show  that  the  agglomerates  are  best  represented  by  sets  of  overlapping  spheres  of 
unequal  radii.  In  this  article,  the  flooded  agglomerate  is  generalized  by  considering  overlapping  spheres 
of  unequal  radii.  As  a  first  cut,  only  two  overlapping  spheres  are  considered.  The  governing  reaction- 
diffusion  equations  are  solved  numerically  using  the  unstructured  finite-volume  method.  The  volumetric 
current  density  is  extracted  for  various  parametric  variations,  and  tabulated.  This  sub-grid-scale  general¬ 
ized  flooded  agglomerate  model  is  first  validated  and  finally  coupled  to  a  computational  fluid  dynamics 
(CFD)  code  for  predicting  the  performance  of  the  PEMFC.  Results  show  that  when  the  agglomerates  are 
small  (<200  nm  equivalent  radius),  the  effect  of  agglomerate  shape  on  the  overall  PEMFC  performance  is 
insignificant.  For  large  agglomerates,  on  the  other  hand,  the  effect  of  agglomerate  shape  was  found  to  be 
critical,  especially  for  high  current  densities  for  which  the  mass  transport  resistance  within  the  agglom¬ 
erate  is  strongly  dependent  on  the  shape  of  the  agglomerate,  and  was  found  to  correlate  well  with  the 
surface-to-volume  ratio  of  the  agglomerate. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

One  of  the  major  impediments  to  commercial  success  of 
hydrogen-oxygen  fuel  cells  is  the  cost  associated  with  the  excessive 
use  of  platinum  within  the  catalyst  layers.  In  a  polymer  electrolyte 
membrane  fuel  cell  (PEMFC),  since  the  catalyst  (platinum)  is  dis¬ 
persed  within  a  complex  porous  matrix  comprised  of  carbon,  the 
ionomer  (Nation)  and  platinum,  the  performance  of  the  fuel  cell 
has  a  convoluted  relationship  with  an  increase  in  the  amount  of 
platinum  within  the  catalyst  layer.  Rather  than  the  amount  of 
platinum,  the  performance  depends  on  how  the  platinum  is  dis¬ 
persed  within  the  porous  matrix  so  that  it  is  effectively  utilized  in 
catalyzing  the  electrochemical  reactions.  Thus,  there  is  strong  moti¬ 
vation  for  optimization  of  the  catalyst  layers  of  a  PEMFC.  Since  the 
proper  functioning  of  a  PEMFC  cathode  requires  existence  of  triple 
phase  boundaries  [1,2]  between  the  ionomer  (for  proton  trans¬ 
fer),  platinum  (for  catalysis)  and  carbon  (for  electron  transfer),  the 
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determination  of  optimum  composition  and  structure  of  the  cata¬ 
lyst  layer  is  a  monumental  task.  Such  studies  are  often  undertaken 
using  experiments  [3,4],  but  are  very  time-consuming  and  expen¬ 
sive.  Computational  modeling  provides  an  alternative  pathway  to 
address  this  critical  issue. 

Model-based  optimization  of  the  cathode  catalyst  layer  of  a 
PEMFC  has  been  a  topic  of  intense  research  over  the  past  two 
decades.  Models  used  for  this  purpose  may  be  broadly  categorized 
as:  (1)  models  that  utilize  the  pseudo-homogeneous  film  concept 
[5-9],  or  the  flooded  agglomerate  concept  [10-28],  and  (2)  direct 
numerical  simulation  (DNS)  of  the  catalyst  layer  [29-33  ].  While  the 
latter  approach  is  more  first-principles  and  general,  it  is  usually 
performed  at  a  much  smaller  scale  and  is  cumbersome  for  cou¬ 
pling  to  a  device-scale  model.  First,  it  requires  reconstruction  of  the 
catalyst  layer  microstructure  from  micrographs.  Secondly,  direct 
numerical  solution  of  the  transport-reaction  equations  within  the 
complex  catalyst  layer  structure  are  difficult  to  perform  and  are 
very  expensive.  Simulation  times  for  a  single  case  may  often  run 
into  days.  This  makes  this  approach  impractical  for  engineering 
calculations.  Nevertheless,  direct  numerical  simulations  are  very 
useful  for  fundamental  understanding  of  the  coupling  between 
transport  and  reactions  at  the  pore  scale  of  the  cathode  catalyst 
layer. 
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Nomenclature 

Af  area  of  face /of  control  volume  (m2) 

Av  total  catalyst  surface  area  per  unit  volume  of  anode 

or  cathode  (nrr 1 ) 

c  dissolved  oxygen  concentration  in  Nation 

(kmolm-3) 

c*  dissolved  oxygen  concentration  in  Nation  in  equi¬ 

librium  with  inlet  gas  (kmol  m-3 ) 
cf)ef  standard  reference  oxygen  concentration 

(kmolm-3) 

Co2ig  oxygen  gas  concentration  in  cathode  gas  pores 
(kmolm-3) 

Cq2  g  oxygen  gas  concentration  at  cathode  inlet 
(kmolm-3) 

d  average  pore  size  of  cathode  (m) 

Vkn  binary  diffusion  coefficient  of  species  k  into  n 

(m2  s-1) 

Dx  diffusion  coefficient  of  water  (m2  s-1 ) 

Dfx  concentration  dependence  of  D^,  dimensionless 

Dj  temperature  dependence  of  Dx,  (m2  s-1 ) 

Do2>n  diffusion  coefficient  of  oxygen  in  Nation  (m2  s-1 ) 
Df2  N  effective  diffusion  coefficient  of  oxygen  in  Nation  in 
agglomerate  (m2  s-1 ) 

F  Faraday’s  constant  (96.487  x  1 06  C  kmol-1 ) 

i  net  current  density  vector  (A  m-2 ) 

iF  ionic  phase  current  density  vector  (A  m-2 ) 

igef  standard  exchange  current  density  on  cathode 

(Am-3) 

is  electronic  phase  current  density  vector  (A  m-2 ) 

i£at  surface  current  density  on  catalyst  surface  of  cath¬ 
ode  (Am-2) 

]k  diffusion  mass  flux  of  the  kth  species  (kg  m-2  s-1 ) 

jjn  net  transfer  current  at  anode  (A  m-3 ) 

j£at  net  transfer  current  density  at  cathode  (A  m-3 ) 

jo  reference  current  density  (A  m-2 ) 

L  cathode  catalyst  layer  thickness  (m) 

mPt  platinum  mass  loading  (kg  m-2 ) 

Mm  molar  mass  of  the  membrane  (kg  kmol-1 ) 

Mk  molecular  weight  of  kth  species  (kg  kmol-1 ) 
n  number  of  electrons  transferred  during  the  electro¬ 

chemical  reaction 

n  number  of  agglomerates  per  unit  volume  of  cathode 

(nr3) 

hf  unit  surface  normal  at  control  volume  face 
N  total  number  of  gas-phase  species 

p  pressure  (Pa) 

Pt|C  platinum-carbon  mass  ratio  in  catalyst  layer  ink, 

dimensionless 

ragg  radius  (actual  or  equivalent)  of  agglomerate  (m) 

r\ ,  r2  radii  of  the  two  overlapping  spheres  (m) 

R  universal  gas  constant  (83 14 J  kmol-1  K-1) 

Sk  production  rate  of  kth  species  (kg  m-3  s-1 ) 

T  absolute  temperature  (I<) 

U  bulk  fluid  velocity  (ms-1 ) 

Vagg  total  volume  of  agglomerate  (m3 ) 

Vnuc  volume  of  nucleus  of  agglomerate  (m3 ) 

Vctg  volume  of  coating  (m3 ) 

Vo  volume  of  cell  or  control  volume  O  (m3 ) 

Yk  mass  fraction  of  kth  species 

Greek 

aa,  ac  Tafel  constants  for  anode,  dimensionless 


aj  Tafel  constant  for  cathode  catalyst  model,  dimen¬ 
sionless 

Pk  concentration  exponents  for  the  kth  species 
8  polymer  coating  thickness  around  agglomerate 

nucleus  (m) 

s  wet  porosity,  dimensionless 

£agg  volume  fraction  of  polymer  in  agglomerate  nucleus, 

dimensionless 

£Cat  porosity  of  cathode  catalyst  layer,  dimensionless 
£S  volume  fraction  of  platinum  +  carbon  in  cathode, 

dimensionless 

£n  volume  fraction  of  polymer  in  cathode,  dimension¬ 
less 

r]  electrode  overpotential  (V) 

pd  electro-osmotic  drag  coefficient,  dimensionless 

k  permeability  (m2) 

A  water  content,  dimensionless 

[Ak]  molar  concentration  of  species  k  (kmolm-3) 

fi  dynamic  viscosity  (kg  m-1  s-1 ) 

p  mass  density  of  mixture  (kg  m-3 ) 

density  of  dry  membrane  (kg  m-3 ) 
a  electrical  conductivity  (£2-1  m-1 ) 

crp  electrical  conductivity  of  the  ionic  phase  ( £2-1  m-1 ) 

gs  electrical  conductivity  of  the  electronic  phase 

(^-1m-1) 

a30  concentration  dependence  of  electrical  conductivity 

(£2-1  m-1) 

0F  ionic  phase  potential  (V) 

0s  electronic  phase  potential  (V) 

4> oc  open  circuit  voltage  (V) 

rcat  tortuosity  of  cathode,  dimensionless 

£  overlap  parameter,  dimensionless 


The  more  popular  approach  is  based  on  hypothesized  mod¬ 
els  of  coupled  mass  transport  and  reactions  within  the  catalyst 
layer  structure.  Historically,  two  different  model  types  have  been 
used  for  this  purpose.  The  first  model  type,  generally  referred  to 
in  the  literature  as  the  pseudo-homogeneous  film  model  [5-9], 
assumes  that  the  catalyst  layer  is  a  porous  matrix  comprised  of 
Nation,  platinum,  and  carbon  in  random  (homogeneous)  config¬ 
uration.  This  model  allows  for  pathways  of  gases,  electrons,  and 
protons  within  the  catalyst  layer,  and  captures  some  of  the  essen¬ 
tial  transport  phenomena  prevalent  in  the  catalyst  layer.  However, 
this  model  does  not  acknowledge  the  necessity  for  the  existence  of 
the  triple-phase  boundary  for  a  functioning  catalyst  layer.  In  con¬ 
trast,  the  flooded  agglomerate  concept,  proposed  in  the  late  1980s 
[10,11],  contends  that  the  platinum  is  supported  on  carbon  parti¬ 
cles,  which  forms  agglomerates  when  mixed  with  an  ionomer.  The 
agglomerate  may  even  be  coated  fully  or  partially  by  an  additional 
ionomer  layer.  The  oxygen  finds  its  way  to  the  platinum  by  first 
dissolving  in  the  ionomer,  and  is  consumed  as  it  transports  to  the 
core  of  the  carbon-platinum  aggregate.  This  model  guarantees  the 
existence  of  triple-phase  boundaries  as  long  as  sufficient  amounts 
of  the  ionomer  are  present.  Calculations  performed  using  this 
model  appear  to  match  experimental  data  better  than  the  pseudo- 
homogeneous  film  model  [8,19].  Most  notable  among  early  studies 
that  have  used  the  flooded  agglomerate  model  is  the  work  of  Jaouen 
etal.  [15,16,34],  Sun  et  al.  [18],andSecanell  et  al.  [19],  although  this 
model  has  found  prolific  usage  more  recently.  Jaouen  et  al.  have  suc¬ 
cessfully  used  this  model  to  predict  the  performance  of  the  cathode 
as  a  function  of  operating  conditions  and  cathode  layer  thicknesses. 
Their  calculations  predicted  experimentally  observed  double 
Tafel  slopes  [16],  attributed  to  local  (within  agglomerate)  mass 
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transport  limitations.  However,  the  model  was  not  exercised  to 
study  effect  of  other  compositional  variables  such  as  Nation  loading 
and  platinum  loading.  Such  studies  were  performed  by  Sun  et  al. 
[18],  who  predicted  an  optimum  Nation  loading— consistent  with 
experimentally  observed  behavior  [3,4].  Studies  performed  by  both 
Sun  et  al.  [  1 8  ]  and  Secanell  et  al.  [  1 9  ]  suggests  that  the  performance 
of  the  cathode  improves  monotonically  as  the  agglomerate  radius  is 
decreased.  However,  it  is  believed  [30,31  ]  that  excessive  reduction 
in  the  agglomerate  radius  can  lead  to  severe  reduction  in  the  size  of 
the  macropores  and  high  tortuosity.  This,  in  turn,  can  result  in  dete¬ 
rioration  of  the  performance,  especially  at  high  current  densities, 
due  to  global  mass  transport  limitations.  A  more  general  form  of 
the  flooded  agglomerate  model  that  accounts  all  of  the  afore-stated 
issues  was  proposed  by  Kamarajugadda  and  Mazumder  [35],  and 
has  been  used  thereafter  by  other  researchers  [27,28]. 

One  of  the  main  assumptions  in  the  flooded  agglomerate  model 
is  that  the  agglomerates  are  spherical  in  shape,  and  that  the  spheres 
are  non-overlapping  or  spatially  isolated.  While  there  is  sufficient 
evidence  from  micrographs  [36]  to  suggest  the  predominance  of 
spherical  shape  of  the  agglomerates  (presumably  because  surface 
tension  tends  to  make  them  spherical),  there  is  little  evidence 
to  suggest  that  these  spheres  are  spatially  isolated.  In  a  recent 
study,  Jain  et  al.  [37]  have  shown  significant  differences  at  the 
agglomerate-scale  in  the  effectiveness  of  oxygen  consumption 
between  plate-like,  spherical,  and  cylindrical  agglomerates.  They 
attribute  this  to  the  differences  in  the  resistance  to  oxygen  diffusion 
between  agglomerates  and  within  the  agglomerates  for  different 
shapes.  Careful  examination  of  micrographs  seem  to  suggest  that 
the  agglomerates  may  be  most  closely  likened  to  groups  of  over¬ 
lapping  spheres  of  unequal  radii,  much  like  clusters  of  soap  bubbles 
that  have  partially  coalesced.  It  is  worth  noting  that  the  assumption 
of  spatially  isolated  spherical  agglomerates  stems  from  the  fact  that 
under  this  assumption,  the  governing  reaction-diffusion  equation 
at  the  agglomerate  scale  can  be  conveniently  solved  using  analytical 
techniques,  which  makes  the  model  easy  to  implement  in  a  large- 
scale  model  for  prediction  of  the  overall  performance  of  a  PEMFC. 
The  assumption  is  not  necessarily  driven  by  observations  of  reality, 
and  therefore,  there  is  a  pressing  need  to  assess  its  validity. 

In  this  study,  the  flooded  agglomerate  model  is  generalized  to 
account  for  agglomerate  shape  effects.  Specifically,  based  on  the 
rationale  presented  earlier  and  as  part  of  this  preliminary  investi¬ 
gation,  agglomerates  represented  by  two  overlapping  spheres  of 
unequal  radii  is  considered,  although  the  model  presented  here 
can  be  used  for  agglomerates  of  arbitrary  shape  and  size.  In  addi¬ 
tion,  two  different  shapes  of  ionomer  coating  are  also  considered. 
The  governing  reaction-diffusion  equation  at  the  agglomerate  scale 
is  solved  numerically  using  the  unstructured  finite-volume  proce¬ 
dure  for  variations  of  relevant  parameters  that  govern  the  catalyst 
layer’s  composition.  The  numerical  solutions  are  post-processed  to 
extract  volumetric  current  densities,  and  are  then  tabulated  for  sub¬ 
sequent  use.  The  agglomerate-scale  model  is  first  validated  against 
analytical  solutions,  and  subsequently  incorporated  in  a  full-scale 


Fig.  2.  Agglomerate  nucleus  represented  by  two  overlapping  spheres  of  unequal 
radii  with  carbon-supported  (black)  catalyst  particles  (red  dots)  distributed  in  the 
polymer  electrolyte  (gray).  (For  interpretation  of  the  references  to  color  in  this  figure 
legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 


computational  fluid  dynamics  (CFD)  code  for  the  prediction  of  the 
performance  of  the  full  PEMFC,  thereby  resulting  in  a  generalized 
multi-scale  model  of  the  PEMFC. 

2.  Model  description 

Computational  modeling  of  a  PEMFC  entails  sub-division  of  the 
entire  fuel  cell  (or  computational  domain)  into  smaller  control  vol¬ 
umes  (or  cells),  and  subsequently  solving  equations  of  conservation 
of  mass,  momentum,  energy,  and  charge  for  each  control  volume. 
The  sub-divided  network  of  control  volumes,  referred  to  as  a  mesh 
or  grid,  is  necessary  to  spatially  resolve  all  quantities  of  interest. 
Since  the  catalyst  layer  of  a  PEMFC  has  features  that  are  less  than 
a  micron  in  size,  while  the  fuel  cell  itself  is  of  the  order  of  a  few 
centimeters  in  size,  it  is  impossible  to  resolve  these  features  within 
the  catalyst  layer  by  the  same  mesh  that  is  used  to  model  the  full 
PEMFC.  Thus,  a  separate  model  is  necessary  at  the  catalyst  layer 
scale,  and  is  referred  to  as  the  sub-grid  scale  model.  This  sub-grid 
scale  model  has  to  be  ultimately  coupled  with  the  larger  scale  model 
in  order  to  study  the  impact  of  small-scale  physics  on  the  device¬ 
scale  performance.  The  relationship  between  the  various  scales  is 
depicted  schematically  in  Fig.  1. 

2.1.  Sub-grid  scale  agglomerate  model 

The  cathode  catalyst  layer  consists  of  three  components: 
the  solid  carbon-platinum  clusters,  the  polymer  electrolyte  (or 
ionomer),  and  gas  pores.  For  the  sub-grid  scale  agglomerate  model 
to  be  valid,  the  agglomerate  size  must  be  sufficiently  small  com¬ 
pared  to  the  cathode  layer  thickness  such  that  the  electric  potential 
within  an  agglomerate  may  be  assumed  to  be  constant.  Trans¬ 
port  of  reactants  from  the  cathode  channel  to  the  cathode  active 
layer  is  dominated  by  diffusion  within  the  pores  of  the  cathode 
diffusion  layer  and  the  cathode  catalyst  layer  (see  Fig.  1).  Oxygen 
reaches  platinum  within  the  agglomerate  by  dissolution  into  the 
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Fig.  1.  Schematic  representation  of  the  cross-section  of  a  PEM  fuel  cell  showing  the  various  length  scales  involved  and  their  relationship. 
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polymer  coating  followed  by  diffusion  through  the  polymer  within 
the  agglomerate.  The  objective  of  the  sub-grid  scale  model  is  to  pro¬ 
vide  the  volumetric  current  density  generated  by  a  unit  volume  of 
the  catalyst  layer  by  taking  into  account  the  effect  platinum  loading, 
Nation  loading,  platinum-to-carbon  ratio,  and  the  microstructure  of 
the  catalyst  layer.  As  discussed  in  Section  1,  the  so-called  flooded 
agglomerate  model  has  been  used  routinely  for  this  purpose. 

In  the  generalized  flooded  agglomerate  model  proposed  here, 
the  cathode  is  assumed  to  be  comprised  of  a  large  number  of 
agglomerates  formed  by  two  overlapping  spheres  of  equal  or 
unequal  radii  (Fig.  2).  In  such  a  scenario,  the  shape  of  the  agglom¬ 
erate  can  be  characterized  by  the  radii  of  the  two  intersecting 
spheres,  r\  and  r2,  and  the  so-called  overlap  parameter,  f.  The 
overlap  parameter  is  another  measure  of  the  distance  between  the 
center  of  the  two  overlapping  spheres,  di2,  which  may  be  written 
as 

di2  =  O  +  £r2  (1) 

where  -1<§<1.§  =  -1  represents  a  single  sphere  of  radius  V\  ( >r2 ), 
while  §  =  1  represents  the  two  spheres  just  touching  each  other. 
Each  agglomerate  is  comprised  of  clusters  of  carbon-supported 
platinum  particles  held  together  by  a  proton-conducting  polymer. 
The  agglomerate  may  be  coated  with  a  thin  film  of  polymer.  In  this 
study,  two  different  coating  shapes  are  considered  (Fig.  3):  (a)  a 
coating  of  uniform  thickness  following  the  contour  of  the  agglom¬ 
erate  nucleus  surface,  and  (b)  a  thin  uniform  coating  around  part 
of  the  curved  nucleus  surface,  with  a  conical  frustum  of  thicker 
layer  connecting  the  two  thin  sections.  Flenceforth,  the  first  coat¬ 
ing  configuration  is  referred  to  as  “thin  coating,”  while  the  second 
configuration  is  referred  to  as  “capsule  coating.”  For  both  these 
coating  configurations,  “coating  thickness  (5)”  refers  to  the  value 
of  thickness  of  the  uniform  coating  around  the  curved  surface  of 
the  agglomerate  nucleus.  Therefore,  for  the  same  prescribed  value 
of  coating  thickness,  the  capsule  coating  case  for  a  given  agglom¬ 
erate  nucleus  has  a  higher  volume  of  polymer  in  the  coating  than 
the  thin  coating  case.  While  the  proposed  model  is  not  in  complete 
agreement  with  observed  microstructures  of  the  catalyst  layer,  it 
represents  the  next  logical  extension  of  an  isolated-sphere  flooded 
agglomerate  model  that  has  been  in  use  thus  far. 

2.1.1.  Agglomerate  model  governing  equation 

Within  the  agglomerate,  oxygen  is  assumed  to  diffuse  through 
the  polymer  electrolyte  to  the  catalyst  surface.  While  it  diffuses 
through  the  nucleus  of  the  agglomerate,  it  reacts  on  the  platinum 
clusters  impregnated  on  the  carbon  surfaces,  and  is  consumed 
(Fig.  2).  For  an  agglomerate  of  arbitrary  shape,  the  reaction- 
diffusion  balance  equation  for  dissolved  oxygen  concentration  in 
Nation,  c,  is  given  by  [  1 0,1 1  ] : 

v  •  (Do^,nvc)  +  fiagg)  =  o  (2) 


where  y  =  1  in  the  nucleus  and  y  =  0  in  the  coating,  Av  represents  the 
total  catalyst  surface  area  per  unit  volume  of  cathode  catalyst  layer 
available  for  the  oxygen  reduction  reaction  (ORR)  [18,19],  i£at  is  the 
surface  current  density  the  catalyst  surface,  and  £a gg  is  the  volume 
fraction  of  polymer  present  within  the  agglomerate  nucleus,  n  is 
the  number  of  electrons  transferred  during  the  ORR  (=2  for  a  H2/02 
PEMFC),  and  F  is  the  Faraday  constant.  The  ORR  kinetics  is  assumed 
to  follow  a  Tafel  law,  and  to  be  first-order  in  oxygen  concentration. 
Hence,  the  surface  current  density  on  the  catalyst  surface  is  given 
by  [15]: 


where  c™*  is  the  standard  reference  oxygen  concentration,  i]jef  is  the 
standard  exchange  current  density  for  a  smooth  catalyst  surface, 
and  r/  is  the  local  overpotential  defined  as  the  difference  between 
the  electronic  and  protonic  phase  potentials,  i.e.,  q  =  </>s  -  0f-  c  and  c* 
are  the  dissolved  oxygen  concentration  in  the  polymer  at  the  cata¬ 
lyst  surface  and  the  dissolved  oxygen  concentration  in  the  polymer 
in  equilibrium  with  the  inlet  gas,  respectively.  Substitution  of  Eq. 
(3)  in  Eq.  (2)  results  in  a  linear  partial  differential  equation  for  the 
dissolved  oxygen  concentration,  c.  N  is  the  effective  oxygen  dif¬ 
fusion  coefficient  in  Nation  within  the  agglomerate,  and  is  obtained 
from  the  Bruggemann  relation  [38]  as: 

D0 9  n  sLg  in  nucleus 

2.  gg  (4) 

Do2,n  in  coating 

The  boundary  condition  for  Eq.  (2 )  must  be  posed  at  the  interface 
with  the  gas  pore.  It  requires  equilibrium  between  the  dissolved 
oxygen  concentration  in  Nation,  c,  and  oxygen  concentration  in  the 
gas  pores,  c02,g,  at  the  outer  surface  of  the  polymer  coating  [15], 
and  is  given  by: 

c*  1 

interface  —  c02,g  =  77c02,g  (5) 

co2,g  n 

where  H  is  Henry’s  constant,  and  cj2  g  represents  the  concentration 
of  oxygen  at  the  cathode  inlet.  At  the  interior  interface  between  the 
coating  and  the  nucleus,  continuity  of  oxygen  concentration  and 
oxygen  flux  must  be  satisfied. 

One  of  the  assumptions  in  the  agglomerate  scale  model  is 
that  the  overpotential  remains  constant  within  the  agglomerate, 
i.e.,  the  ionic  and  electronic  phase  potentials  do  not  vary  within 
the  agglomerate.  This  may  be  justified  using  a  simple  order-of- 
magnitude  analysis.  The  potential  drop  from  the  center  to  the 
surface  of  a  spherical  agglomerate  of  diameter  1 000  nm  (the  largest 
one  considered  here)  using  Ohm’s  Law  is  given  by  Arp^ir-i/a, 
where  i  is  the  current  density,  r\  is  the  radius  of  the  agglomerate, 
and  a  is  either  the  electronic  or  protonic  conductivity,  depending 
on  whether  the  protonic  phase  potential  drop  or  the  electronic 
phase  potential  drop  is  being  considered.  Using  typical  values  of 
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i  =  104Airr2,  r\  =  500 nm,  and  <r  =  10Sm1  (for  Nafion),  we  get 
A </>  =  0.05  mV.  This  is  the  maximum  possible  potential  drop  in  the 
ionic  phase  potential  under  any  condition.  For  electron  transport 
the  potential  drop  is  one  order  of  magnitude  lower.  Since  0.05  mV 
is  negligible  (4  orders  of  magnitude  lower)  compared  to  the  cell 
voltage  of  ~1  V,  it  can  be  safely  neglected. 


2.1.2.  Numerical  solution  of  agglomerate  model  equation 

In  order  to  solve  Eq.  (2)  for  an  agglomerate  of  arbitrary  shape, 
the  finite-volume  method  is  employed.  In  this  method,  the  entire 
agglomerate  is  first  discretized  into  small  control-volumes.  In  this 
work,  an  unstructured  mesh  is  used  to  allow  flexibility  of  modeling 
any  arbitrary  agglomerate  shape.  The  unstructured  mesh  is  gen¬ 
erated  using  the  commercial  mesh  generation  code  CFD-GEOM™. 
In  the  finite-volume  method,  the  governing  equation  [Eq.  (2)]  is 
first  integrated  over  a  representative  control  volume,  O.  This  step 
is  followed  by  application  of  the  Gauss-divergence  theorem,  which, 
following  procedures  outlined  elsewhere  [39,40],  yields  a  set  of 
discrete  algebraic  equations  of  the  following  form 


eff 

o2.n  / 


(  CN,/  -  C0 


.  h/Af 


Vo  =  S0V o 


(6) 


where  the  summation  is  over  all  faces  of  the  control  volume,  O.  V0 
is  the  volume  of  the  control  volume  in  question,  and  Af  is  the  sur¬ 
face  area  of  each  face  of  the  control  volume.  D? ff  M ,  is  the  diffusion 
coefficient  of  oxygen  at  the  face,  and  is  obtained  from  cell-center 
diffusion  coefficients  using  distance-weighted  interpolation.  cN^is 
the  concentration  at  the  cell  center  on  the  other  side  of  face/,  i.e., 
the  neighboring  cell  N,  c0  is  the  concentration  at  cell  O,  and  fy  is 
the  distance  between  the  cell  centers  of  the  two  cells  straddling 
face /in  the  direction  of  the  surface  normal.  JTj  represents  the  flux 
tangential  on  the  face/,  and  the  procedure  to  express  it  in  terms  of 
the  cell-center  concentrations  may  be  found  elsewhere  [39,40]. 

Eq.  (6)  represents  a  set  of  coupled  linear  algebraic  equations 
whose  solution  yields  the  cell-center  concentrations  at  all  cells 
within  the  entire  agglomerate.  This  system  of  coupled  algebraic 
equations  was  solved  using  the  Generalized  Minimum  Residual 
Solver  (GMRES)  [41]  with  incomplete  LU  (ILU)  pre-conditioning. 
Typically,  20  Krylov  sub-spaces  were  used,  and  the  residual 
was  reduced  by  six  orders  of  magnitude.  Example  of  a  sample 
unstructured  mesh  and  the  computed  solution  (dissolved  oxygen 
concentration)  is  depicted  in  Fig.  4.  This  particular  calculation  was 
performed  for  an  agglomerate  of  a  pair  of  intersecting  spheres  with 
r\  =1000  nm,  r2  =  600  nm,  an  overlap  parameter  of  f  =  0,  and  with 
a  coating  of  thickness  5  =  60  nm.  It  is  clear  from  the  figure  that  the 
mesh  resolution  used  for  such  computations  is  quite  high  so  as 
to  accurately  capture  the  oxygen  concentration  distribution,  espe¬ 
cially  in  the  region  where  the  spheres  intersect.  Grid  independence 
studies  revealed  that  any  mesh  with  less  than  500,000  cells  is  inad¬ 
equate  for  such  computations.  In  terms  of  computational  effort,  the 
mesh  generation  process  was  the  primary  consumer  of  CPU  time. 
The  mesh  generation  for  each  mesh  required  about  45  minutes  of 
CPU  time  on  a  3.2  GEIz  Intel  Pentium4  processor,  while  the  actual 
solution  of  the  governing  equation  required  less  than  10  minutes 
in  most  cases. 


Fig.  4.  Results  from  a  demonstration  case  of  the  numerical  solution  of  the  gen¬ 
eralized  flooded  agglomerate  model  equation:  (a)  a  typical  unstructured  mesh 
comprised  of  784,029  cells  and  (b)  contours  of  02  concentration  distribution  (in 
kmol  irr3 ). 


In  order  to  compute  the  effective  properties  for  transport  in  the 
cathode  catalyst  layer,  the  volume  fractions  of  the  carbon-platinum 
solid  clusters  ($5).  the  proton-conducting  polymer  (eN),  and  the 
pores  (£Cat),  need  to  be  calculated.  Also,  the  total  catalyst  surface 
area  per  unit  volume  of  cathode  catalyst  layer  available  for  the 
oxygen  reduction  reaction  (ORR),  Av,  needs  to  be  computed.  These 
quantities  are  computed  using  the  relationships  presented  in  ear¬ 
lier  publications  [18,19,35],  and  are  omitted  here  for  the  sake  of 
brevity.  Under  the  assumption  that  the  agglomerate  nucleus  is 
made  up  of  only  the  solid  component  and  the  polymer,  the  num¬ 
ber  of  agglomerates  per  unit  volume  required  to  produce  the  solid 
component  volume  fraction,  es»  is  given  by: 


Wiuc(l  —  £agg) 


2.1.3.  Computation  of  volume -average  quantities 

Once  the  solution  to  the  agglomerate-scale  governing  equation 
[Eq.  (6)]  has  been  obtained,  it  can  be  post-processed  to  extract 
quantities  that  appear  as  inputs  to  the  device-scale  CFD  model,  to 
be  described  in  the  next  section. 


where  Vnuc  is  the  volume  of  the  agglomerate  nucleus,  and  is  com¬ 
puted  directly  by  the  agglomerate  equation  solver  by  summing  up 
the  volumes  of  all  the  control  volumes  that  constitute  the  nucleus. 
Further,  assuming  that  the  entire  polymer  electrolyte  (Nafion)  in 
the  cathode  catalyst  layer  is  only  present  either  in  the  agglomer¬ 
ate  nucleus  or  in  the  polymer  coating  around  the  agglomerate,  the 
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volume  fraction  of  the  polymer  electrolyte  in  the  cathode  can  be 
obtained  using 

£n  =  n(V nuc^agg  +  ^ctg)  (8) 

where  Vctg  is  the  volume  of  the  polymer  coating  around  the  agglom¬ 
erate  nucleus,  and  is  also  computed  directly  by  the  agglomerate 
equation  solver.  Once  the  volume  fractions  of  the  solid  component 
and  the  polymer  are  obtained,  the  porosity  of  the  cathode  catalyst 
layer  is  obtained  using 


Scat  —  1  —  £s  — 


(9) 


The  current  per  unit  volume  generated  by  the  agglomerate,  j^gg, 
is  calculated  using: 


;agg 
J  T 


Ivaee  “  £agg)dV 


1  In 


= 


i-ref 


^agg 


rref 
°2  > 


x  exp 


A/(l  -  £agg)  fv nuc  CC ^ 


^agg 


(10) 


where  Vagg  is  the  total  volume  of  the  agglomerate  (nucleus  plus 
coating).  The  integral  in  the  numerator  of  Eq.  (10)  is  computed 
directly  by  the  agglomerate  equation  solver  by  summing  over  the 
volume  weighted  concentration  of  each  control  volume.  The  cur¬ 
rent  generated  per  unit  volume  of  the  cathode  may  be  expressed  in 
terms  of  the  current  generated  per  unit  volume  of  the  agglomerate 
as 


jcat  =  (1_£cat)j-agg  (n) 

As  stated  earlier,  the  quantity  j£at  appears  as  a  source  in  the 
governing  equations  for  current  conservation  within  the  cathode, 
as  will  be  discussed  shortly. 


term  due  to  electrochemical  reactions,  and  is  non-zero  only  in  the 
catalyst  layers.  It  results  from  the  fact  that  a  conservation  equation 
for  protons  is  not  directly  solved,  and  therefore,  the  mass  of  protons 
created  or  destroyed  need  to  be  added  or  subtracted  out  of  the 
overall  continuity  equation.  The  last  term  in  Eq.  (13)  represents  the 
sub-grid  scale  drag  force  imposed  by  the  pore  walls  on  the  fluid, 
written  using  the  linear  Darcy’s  Law  [38].  In  purely  open  regions, 
such  as  in  the  gas  channels,  e  ->  1  and  k  ->  oo,  and  Eqs.  (12)  and  (13) 
reduce  to  the  well-known  Navier-Stokes  equations.  In  Eq.  (14),  Yk 
is  the  mass  fraction  of  the  kth  species,  J/<  is  the  mass  diffusion  flux 
of  the  kth  species,  and  Sk  is  the  production  rate  of  the  kth  species 
due  to  electrochemical  reactions.  The  total  number  of  species  in  the 
system  is  denoted  by  N. 

The  diffusion  flux  of  the  kth  species,  ]k,  is  modeled  using  the 
so-called  dilute  approximation  [46].  By  this  approximation,  the 
diffusion  flux  is  written  as 


J  k  —  PDkm^Yk 


(15) 


where  Dkm  is  the  free-stream  diffusivity  of  species  k  into  the 
mixture,  and  is  henceforth  denoted  by  Dk  for  simplicity.  The  free- 
stream  diffusivity  is  given  by  the  relation  [46]: 


Dkm  ~  Dk  — 


1-Xk 


i  #  k 


(16) 


where  Vki  is  the  binary  diffusivity  of  species  k  into  species  i,  and  Xk 
is  the  mole  fraction  of  the  kth  species.  Substitution  of  Eq.  (15)  into 
Eq.  ( 1 4)  yields  the  appropriate  species  transport  equation  under  the 
dilute  approximation  formulation  for  open  regions: 

V  •  ( PUYk )  =  V  •  ( pDkVYk )  +  Sk  Vk  =  1 ,  2, . . . ,  N  (17) 


2.2.  Computational  fluid  dynamics  (CFD)  model 

The  sub-grid  scale  catalyst  layer  model,  described  in  the 
preceding  section,  was  implemented  into  a  two-dimensional  com¬ 
putational  fluid  dynamics  (CFD)  code  that  has  been  developed 
specifically  for  PEMFC  applications,  and  has  been  validated  [42] 
and  successfully  used  in  previous  studies  [35].  The  purpose  of  this 
coupling  is  to  enable  prediction  of  the  overall  performance  (polar¬ 
ization  behavior)  of  the  fuel  cell  by  taking  into  account  both  local 
(agglomerate-scale)  and  global  (device-scale)  losses.  The  governing 
equations  are  the  equations  of  conservation  of  mass  (both  overall 
and  individual  species),  momentum,  energy,  and  charge.  For  sim¬ 
plicity,  it  is  assumed  that  the  temperature  does  not  change  spatially, 
and  thus,  the  energy  conservation  equation  is  not  solved.  Further¬ 
more,  it  is  assumed  that  water  exists  only  in  its  vapor  phase  in  gas 
diffusion  layers  and  channels,  and  two-phase  effects  are  not  con¬ 
sidered.  The  governing  conservation  equations  are  different  inside 
and  outside  the  membrane  and  are,  therefore,  presented  here  sep¬ 
arately. 

2.2.1.  Channels,  gas  diffusion  layers  (GDL),  and  active  catalyst 
layers  (ACL) 

The  governing  conservation  equations  for  mass  and  momentum, 
at  steady  state,  are  written  as  [43-45]: 

overall  mass:  V  •  (spU)  =  5m  (12) 

momentum  :  V  •  (peUU)  =  -eVp  +  V  •  (pteffVU)  +  (13) 

species  mass  :  V  ■  (pUYk)  = -V  Jk +  Sk  Vk  =  l,2,...,N  (14) 

where  p  is  the  density,  p  is  the  pressure,  /z  is  the  dynamic  viscosity 
of  the  fluid,  and  U  is  the  fluid  velocity  vector.  £  is  the  porosity  of  the 
medium  and  /c  is  the  permeability.  In  Eq.  (12),  Sm  is  the  mass  source 


In  porous  regions,  it  is  customary  to  use  the  so-called 
Bruggemann  relation  [38],  and  modify  the  free-stream  diffusion 
coefficient,  such  that  the  governing  species  conservation  equation 
becomes 


V  •  (£pUYk)  =  V  •  ( pDksTVYk )  +  Sk  Vk  =  1, 2, . . . ,  N  (18) 


where  r  is  the  tortuosity  of  the  porous  region.  Traditionally,  a  value 
of  r  =  1 .5  is  used  in  Eq.  ( 1 8).  In  order  to  examine  the  effect  of  cathode 
structure  and  tortuosity  on  fuel  cell  performance,  a  more  realistic 
tortuosity  model,  proposed  by  Abbasi  et  al.  [47],  is  used.  Hence  the 
governing  species  conservation  equation  in  the  cathode  catalyst 
layer  becomes: 

v  •  (ecatpWfc)  =  V  .  (pDfcgvYk)  +Sk  Vk  =  l,2,...,N  (19) 

where  the  tortousity,  rcat*  is  given  by  [47] 

teat  = —  +  1.196^  (20) 

Scat  d 

where  crdev  is  the  standard  deviation  of  the  pore  size  in  the  cathode 
catalyst  layer.  In  this  study,  a  constant  value  of  crdev  =  50  nm  has 
been  used  for  the  sake  of  convenience,  d  is  the  average  pore  size  of 
the  cathode  catalyst  layer,  and  is  given  by  [48]: 


d 


4  Scat 

3(1^0  agS 


(21) 


where  ragg  is  the  radius  of  the  equivalent  sphere  that  has  the  same 
volume  as  the  overlapping  spheres  used  here. 

The  source  due  to  the  electrochemical  reactions  is  non-zero  only 
in  the  active  catalyst  layers  of  the  anode  and  cathode,  and  is  zero 
elsewhere.  It  is  written  as  [49] 


Sk  = 


jfMk 
2  F 


(22) 
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where  jj  is  the  net  transfer  current  due  to  electrochemical  reac¬ 
tion,  Mk  is  the  molecular  weight  of  the  kth  species,  and  F  is  the 
Faraday  constant.  The  anode  transfer  current  jf1  is  expressed  using 
the  global  single-scale  Butler-Volmer  equation  [49]: 


although  there  is  some  evidence  that  pressure-driven  advection 
may  also  occur  [50-54].  In  phenomenological  membrane  models, 
it  is  customary  to  express  water  transport  in  terms  of  the  water 
content,  X,  as  [55]: 


Jt"  =joA™ 


exp(¥^)-exp(-¥^)]n[^iA 

k= 1 


(23) 


where  j0  is  the  reference  current  density,  aa  and  ac  are  kinetic 
constants  determined  from  experimentally  measured  Tafel  plots. 
[Ak]  and  fik  are  the  molar  concentrations  and  the  concentra¬ 
tion  exponents  for  the  kth  species,  respectively,  ij  is  the  electrode 
overpotential  (including  both  activation  and  concentration  overpo¬ 
tential)  and  is  defined  as  the  difference  between  the  electronic  or 
solid  phase  potential  (</>s)  and  the  electrolyte  or  ionic  phase  poten¬ 
tial  (0F),  z.e.,  tj  —  0s  —  0F*  In  this  study,  the  above  Butler-Volmer 
kinetic  relation  was  used  only  for  the  anode  and  no  multi-scale 
model  was  used,  while  the  cathode  transfer  current  was  obtained 
from  Eq.  (11),  following  the  multi-scale  procedure  outlined  in  Sec¬ 
tion  2.1. 

In  a  PEMFC  calculation,  in  addition  to  mass  conservation,  it  is 
necessary  to  enforce  charge  conservation.  Under  the  assumption 
of  electro-neutrality,  charge  conservation  reduces  to  current  con¬ 
servation,  written  as 


V  •  i  =  0  (24) 

where  i  is  the  current  density  vector.  In  a  PEMFC,  the  current  flow  is 
due  to  protons  (H+)  flowing  through  the  membrane,  resulting  in  an 
ionic  phase  current  (iF),  and  due  to  electrons  flowing  through  the 
carbon  in  the  porous  matrix  of  the  gas  diffusion  layers,  resulting  in 
an  electronic  phase  current  (is).  Thus,  Eq.  (24)  can  be  rewritten  as 
follows: 


V  •  iF  +  V  •  is  =  0  (25) 

Since  current  transport  in  the  ionic  phase  is  due  to  ions  and  that 
in  the  electronic  phase  is  due  to  electrons,  the  transport  in  each 
phase  is  governed  by  separate  electric  potential  fields.  Using  Ohm’s 
law,  Eq.  (25)  can  be  written  as 

V  •  (crFV0F)  +  V  •  (<7sV0s)  —  0  (26) 

where  aF  and  <rs  are  the  conductivities  of  the  ionic  and  electronic 
phases,  respectively.  The  exchange  of  current  from  the  ionic  to  the 
electronic  phase  occurs  due  to  electrochemical  reactions  during 
which  electrons  are  transferred  from  one  phase  to  the  other.  Thus, 
Eq.  (26)  can  be  rewritten  as  [44,49]: 

-V  •  (crFV0F)  =  V  •  (crsV0s)  =Jt  (27) 

Eq.  (27)  represents  two  elliptic  partial  differential  equations  that 
are  strongly  coupled  through  the  transfer  current  source.  The  ionic 
phase  electric  potential  equation  (for  0F)  must  be  solved  in  the 
active  catalyst  layer  and  the  membrane,  while  the  electronic  phase 
electric  potential  equation  (for  0S)  must  be  solved  in  the  active 
catalyst  layer  (ACL)  and  the  gas  diffusion  layer  (see  Fig.  1).  In  the 
active  catalyst  layer,  both  equations  are  solved,  and  are  strongly 
coupled.  The  difference  in  value  between  0S  and  0F  represents  the 
total  electrode  overpotential. 


2.2.2.  Water  and  current  transport  in  membrane 

Nation  membranes  are,  for  all  practical  purposes,  impermeable 
to  all  gases  except  water.  Thus,  the  governing  equation  for  species 
transport  outside  the  membrane  [Eq.  (14)]  is  irrelevant  in  this  con¬ 
text.  The  only  species  that  needs  to  be  considered  is  water.  For  all 
other  species,  zero  flux  boundary  conditions  must  be  used  at  the 
membrane-ACL  interface.  Water  transport  in  the  membrane  of  a 
PEMFC  occurs  primarily  due  to  diffusion  and  electro-osmotic  drag, 
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Pm 

Mm 


D\  VA 


(28) 


The  water  content  is  defined  as  the  ratio  of  the  number  of  water 
molecules  to  the  number  of  SO3-  charge  sites  in  the  medium.  In 
Eq.  (28),  it  has  been  assumed  that  advective  transport  is  negligible. 
Dx  denotes  the  diffusion  coefficient  of  water  in  Nation  expressed 
in  terms  of  the  water  content,  pd  denotes  the  electro-osmotic  drag 
coefficient,  and  pm/Mm  represents  the  molar  density  of  the  mem¬ 
brane.  Fuller  and  Newman  [56]  give  the  diffusion  coefficient  as: 


Dx  =  (2.1  x  10  7)  x  X  x  exp 


’  2436' 
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(29) 


where  T is  the  absolute  temperature.  The  electro-osmotic  drag  coef¬ 
ficient  is  given  by  [57]: 

riA  =  (30) 

Current  transport  in  the  membrane  is  described  by  the  ionic 
phase  part  of  Eq.  (27),  except  that  there  is  no  source,  z.e.,  the  right 
hand  side  of  the  equation  is  zero.  The  ionic  phase  electrical  con¬ 
ductivity  of  the  membrane,  <rF,  is  again  expressed  by  the  empirical 
correlation  [57]: 


crF  =  100  exp 


1268  (  JL 

\  303 


(31) 


where  cr30  is  the  electrical  conductivity  of  the  membrane  at  30  °C, 
and  is  given  by  [56]: 


<x3o  =  0.0051391 -0.00326  fori  >  1 


(32) 


The  ionic  phase  electrical  conductivities  in  the  anode  cata¬ 
lyst  layer  and  the  cathode  catalyst  layer  are  evaluated  using  Eqs. 
(31)  and  (32)  based  on  local  values  of  water  content.  However, 
in  these  regions,  the  conductivities  need  to  be  corrected  in  order 
to  account  for  the  fact  that  the  proton-conducting  polymer  is  just 
one  of  the  components  in  the  layer.  The  effective  ionic  phase  elec¬ 
trical  conductivity  in  the  cathode  is  given  by  Jaouen  et  al.  [15], 
wherein  non-overlapping  spherical  agglomerates  are  assumed.  For 
the  present  study,  the  expression  provided  by  Jaouen  et  al.  [15]  is 
modified  to  yield 
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where  the  constant  a\  has  been  introduced  following  Mazumder 
and  Kamarajugadda  [35]  to  ensure  that  in  the  case  when  the 
agglomerates  do  not  have  any  coating,  the  effective  protonic  con¬ 
ductivity  decays  to  zero.  The  constant  a\  is  given  by 


a\  =  min 


0, 


Wiuc  / 


1/3 


+  (1 


(34) 


It  is  evident  from  Eqs.  (28)-(32)  that  the  two  quantities  that 
dictate  the  mass  transport  and  electrical  properties  of  a  Nation 
membrane  are  its  water  content  and  temperature.  The  water  con¬ 
tent  is  obtained  by  solution  of  the  conservation  equation  for  water 
[Eq.  (28)].  However,  solution  of  this  equation  requires  specifica¬ 
tion  of  either  the  water  content  itself  or  the  flux  of  water  at  the 
membrane-ACL  interface,  which  in  turn  requires  coupling  with 
the  rest  of  the  computational  domain.  Various  approaches  for  cou¬ 
pling  the  membrane  with  the  overall  calculation  procedure  and  for 
implementation  of  the  interface  conditions  are  discussed  in  Kama¬ 
rajugadda  and  Mazumder  [42]. 
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Fig.  5.  Variation  with  agglomerate  radius  and  polymer  coating  thickness  of  (a) 
numerically  computed  current  generation  per  unit  volume  (in  Am-3)  by  spheri¬ 
cal  agglomerates  and  (b)  percentage  difference  between  numerical  and  analytical 
solution. 

3.  Results  and  discussion 

Prior  to  performing  full-scale  CFD  calculations  to  predict  the  per¬ 
formance  of  a  PEMFC,  validation  of  the  numerical  solution  of  the 
agglomerate  scale  equation  was  deemed  necessary.  Following  this 
important  step,  the  sub-grid  scale  model  was  implemented  into 
the  large-scale  CFD  model,  and  the  effect  of  the  agglomerate  shape 
on  the  overall  performance  of  the  fuel  cell  was  investigated.  These 
studies  are  presented  next. 

3.1.  Validation  of  the  agglomerate  scale  model 

The  purpose  of  the  sub-grid  scale  agglomerate  model  is  to  pro¬ 
vide  a  measure  of  the  volumetric  current  density  generated  by  the 
agglomerates.  Therefore,  the  numerical  solution  of  the  generalized 
flooded  agglomerate  model  equation  is  validated  against  the  ana¬ 
lytical  solution  [15]  for  the  current  density  generated  by  spherical 
agglomerates,  before  using  it  to  analyze  the  effect  of  non-spherical 
agglomerate  size  and  shape  on  overall  fuel  cell  performance. 
Fig.  5(a)  shows  the  volumetric  current  generation  predicted  by 
the  numerical  solution  of  the  generalized  flooded  agglomerate 
model  equation  [Eq.  (2)]  for  a  large  range  of  agglomerate  radius 
and  polymer  coating  thickness  for  nominal  cathode  overpoten¬ 
tial  value  of  r)  =  -0.5  V,  platinum  loading  of  mpt  =  0.45  mg  cm-2, 
Nation  volume  fraction  of  £agg  =  0.4,  and  platinum-carbon  ratio 
of  Pt|C  =  0.28.  Fig.  5(b)  depicts  the  percentage  difference  between 
the  numerical  solution  and  the  analytical  solution  [15].  It  is  clear 
from  Fig.  5(b)  that  the  numerical  solutions  are  within  ±1%  of  the 


Table  1 

Values  of  cathode  model  constants  and  properties. 


Model  parameter 

Value/method  of  calculation 

Transfer  coefficient  at  cathode 

1  (>0.8  V)  [18] 

(Tafel  constant) 

0.55  (<0.8  V)  [18] 

Reference  exchange  current 

1.288  x  10-4  Acm-2  (>0.8  V)  [18] 

density  at  cathode 

2.291  x  10-2  Acm-2  (<0.8  V)  [18] 

Diffusion  coefficient  of  O2  in  Nation 

3.966  x  10-10m2s-1  [18] 

Henry’s  constant 

10  [18] 

analytical  solutions  for  the  entire  range  of  agglomerate  radius  and 
coating  thickness  considered  in  this  study.  Although,  the  above 
results  do  not  automatically  prove  the  reliability  of  the  numeri¬ 
cal  procedure  for  non-spherical  shapes,  based  on  the  results  of  this 
single  validation  study,  the  generalized  flooded  agglomerate  model 
can  be  used  with  some  amount  of  confidence  to  predict  current 
density  generated  by  agglomerates  of  arbitrary  shapes. 

3.2.  Effect  of  agglomerate  shape  and  size  on  local  current 
generation 

Following  the  validation  studies,  the  generalized  flooded 
agglomerate  model  was  used  to  clearly  elucidate  the  effect  of 
agglomerate  shape  and  size  on  the  current  density  generated 
by  agglomerates  characterized  as  pairs  of  intersecting  spheres. 
Agglomerate  size  was  varied  using  three  different  values  (200  nm, 
600  nm,  and  lOOOnm)  for  the  two  radii  (rlf  r2)  of  the  spheres. 
This  gives  rise  to  six  different  agglomerate  configurations— three  of 
which  have  equal  sized  overlapping  spheres,  and  three  of  unequal 
sized  overlapping  spheres.  For  each  of  the  above  six  cases,  the 
agglomerate  shape  is  further  influenced  by  the  distance  between 
the  two  centers,  which  can  be  characterized  by  the  overlap  param¬ 
eter  (§).  In  this  study,  the  overlap  parameter  was  varied  from  -1  to 
+1 .  Also,  in  order  to  examine  the  effect  of  polymer  coating,  calcula¬ 
tions  for  the  current  generated  per  unit  volume  for  each  of  the  above 
six  configurations  are  conducted  for  agglomerates  with  no  coating, 
with  thin  coating,  and  with  capsule  coating.  The  coating  thickness 
for  each  agglomerate  was  calculated  based  on  the  radius  of  the 
larger  of  the  two  intersecting  spheres  using  a  value  of  8/r^  =  0.06. 
All  calculations  discussed  in  this  sub-section  are  carried  out  with  a 
nominal  cathode  overpotential  value  of  p  =  -0.5  V.  The  thickness  of 
the  cathode  is  assumed  to  be  L  =  15  pan,  the  polymer  volume  frac¬ 
tion  within  the  agglomerate  nucleus  is  e  =  0.4,  the  platinum  mass 
loading  is  mPt  =  0.45  mg  cm-2,  and  the  platinum-to-carbon  ratio  in 
the  catalyst  ink  is  Pt|C  =  0.28.  Other  relevant  parameters  used  for 
these  calculations  are  listed  in  Table  1.  All  calculations  were  per¬ 
formed  assuming  that  the  oxygen  concentration  in  the  gas  pore 
surrounding  the  agglomerate  to  be  equal  to  the  oxygen  concen¬ 
tration  at  the  fuel  cell  cathode  flow  channel  inlet.  Evaluation  of 
current  generated  per  unit  volume  of  agglomerate  for  other  values 
of  oxygen  concentration  is  straightforward  because  of  the  linear 
dependence  of  the  current  density  on  the  concentration  distribu¬ 
tion,  as  is  evident  from  Eq.  (10). 

Since  the  number  of  parametric  studies  conducted  [58]  is  too 
large  to  discuss  extensively,  only  the  important  finding  from  this 
study  is  discussed  here.  The  studies  revealed  that  agglomerates 
of  intersecting  spheres  resulted  in  a  higher  current  generated  per 
unit  volume  in  comparison  to  spherical  agglomerates  of  equiva¬ 
lent  volume.  The  current  generated  per  unit  volume  was  found  to 
correlate  well  with  the  surface-area-to-volume  ratio  of  the  agglom¬ 
erate,  and,  for  a  given  volume,  intersecting  spheres  have  a  higher 
surface-area-to-volume  ratio  than  spherical  agglomerates.  This  is 
depicted  in  Fig.  6.  In  this  figure,  the  normalization  is  with  respect 
to  a  single  sphere  (the  larger  of  the  two  radii,  i.e.,  r\ ).  It  is  clear  from 
Fig.  6(b),  that  the  current  generation  improves  for  an  agglomerate 
comprised  of  intersecting  spheres  in  comparison  to  a  single  sphere 
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Fig.  6.  Effect  of  agglomerate  size  and  shape  on  (a)  surface-area-to-volume  ratio  and 
(b)  local  current  generated  per  unit  volume. 


of  the  same  volume  (“equivalent  sphere”).  The  improvement  is 
more  pronounced  for  larger  spheres  than  smaller  ones  because  in 
smaller  spheres,  the  mass  transport  resistances  are  not  significant 
to  begin  with,  and  the  entire  agglomerate  is  almost  saturated  with 
the  reactant  (02).  Further,  comparison  of  Fig.  6(a)  with  Fig.  6(b) 
clearly  illustrates  that  the  current  density  correlates  strongly  with 
the  surface-area-to-volume  ratio  of  the  agglomerate,  although  the 
relationship  is  not  quite  linear.  While  the  results  presented  in  Fig.  6 
are  for  a  case  without  coating,  similar  trends  were  obtained  with 
both  capsule  and  thin  coating. 

3.3.  Effect  of  agglomerate  shape  and  size  on  overall  PEMFC 
performance 

In  addition  to  affecting  the  current  density  generated  at  the 
agglomerate  scale,  agglomerate  shape  and  size  affect  global  prop¬ 
erties  like  cathode  porosity,  tortuosity,  and  proton  conductivity, 
thereby  affecting  the  overall  fuel  cell  performance  in  more  ways 
than  one.  For  example,  one  of  the  findings  of  the  studies  described 
in  the  previous  section  is  that  the  current  generated  by  agglomer¬ 
ates  without  coating  is  larger  than  that  generated  by  agglomerates 
with  coating  because  of  reduced  mass  transport  resistance.  How¬ 
ever,  lack  of  coating  can  cause  difficulty  in  proton  transport  within 
the  catalyst  layer  due  to  poor  contact,  as  discussed  in  Section  2.2.2. 
Therefore,  the  results  of  the  agglomerate  scale  modeling  cannot  be 


Table  2 

Look-up  table  with  current  per  unit  volume  generated  by  agglomerates  of 
200-200  nm  intersecting  spheres  with  £  =  0  for  various  values  of  cathode 
overpotential. 


Overpotential  (V) 

Volumetric  current  density  (Am-3) 

No  coating 

Thin  coating 

Capsule  coating 

0 

4.34852  xlO2 

3.72223  xlO2 

3.60637  x  102 

-0.05 

2.62150  x  103 

2.24394  x  103 

2.17409  x  103 

-0.1 

1.58036  xlO4 

1.35275  xlO4 

1.31065  xlO4 

-0.1404 

6.74781  x  104 

5.77596  x  104 

5.59617  x  104 

-0.1904 

4.06774  x  105 

3.48187  x  105 

3.37349  x  105 

-0.2404 

2.45160  x  106 

2.09847  x  106 

2.03313  x  105 

-0.2904 

2.37045  x  107 

2.02860  xlO7 

1.96534  x  107 

-0.3404 

6.31206  xlO7 

5.39983  xlO7 

5.23088  x  107 

-0.3904 

1.66735  x  108 

1.42502  x  108 

1.38006  x10s 

-0.4404 

4.31615  xlO8 

3.67993  xlO8 

3.56142  x  108 

-0.4904 

1.06647  xlO9 

9.04054  xlO8 

8.73617  x  108 

-0.6 

5.50848  x  109 

4.49353  x  109 

4.31027  x  109 

-0.7 

1.72047  x  1010 

1.27908  xlO10 

1.21477  x  1010 

-0.8 

4.58003  x  1010 

2.80301  xlO10 

2.63290  x  1010 

-0.9 

9.71473  x  1010 

4.52061  x  1010 

4.21920  x  1010 

-1 

1.44451  x  1011 

5.51051  x  1010 

5.13459  x  1010 

directly  extrapolated  to  overall  PEMFC  performance.  To  gain  a  bet¬ 
ter  understanding  of  the  overall  PEMFC  performance,  the  results 
of  the  agglomerate  scale  model  must  be  used  in  conjunction  with 
larger  length  scale  effects. 

Unlike  the  spherical  flooded  agglomerate  model  where  a  closed- 
form  analytical  expression  for  the  volumetric  current  generated 
by  an  agglomerate  is  available,  the  generalized  flooded  agglom¬ 
erate  model  requires  calculation  of  the  current  density  based  on 
the  numerical  solution  of  the  governing  equation  [Eq.  (2)].  In  other 
words,  the  generalized  flooded  agglomerate  model  provides  data 
for  current  generated  per  unit  volume  of  agglomerate  for  discrete 
combinations  of  the  cathode  overpotential  (77),  and  the  oxygen  con¬ 
centration  ( Cq2  g )  in  the  gas  pore  surrounding  the  agglomerates, 
even  with  the  other  cathode  model  parameters  being  held  con¬ 
stant.  Therefore,  scaling  and  interpolation  are  necessary  in  order 
to  obtain  current  density  values  for  intermediate  values  of  r\  and 
Cq2  g  that  are  encountered  during  iterations  in  the  numerical  solu¬ 
tion  of  the  governing  equations  within  the  CFD  model  for  PEMFCs. 
Accordingly,  for  each  of  the  agglomerate  configurations  listed  at  the 
beginning  of  this  sub-section,  the  generalized  flooded  agglomerate 
model  was  used  to  first  create  a  look-up  table  of  values  of  current 
generated  per  unit  volume  [using  Eq.  (10)]  for  various  specified 
values  of  r\  ranging  from  zero  to  the  open  circuit  potential,  with 
an  oxygen  concentration  equal  to  the  concentration  of  oxygen  at 
the  fuel  cell  cathode  flow  channel  inlet  (cj2  g).  Since  the  equation 
governing  the  oxygen  reduction  reaction  within  the  agglomerate 
[Eq.  (2)]  is  linear,  the  cathode  current  density  at  an  arbitrary  con¬ 
centration  of  oxygen  in  the  gas  pore  (c02ig)  can  be  simply  obtained 
by  scaling.  Since  the  overpotential  appears  in  the  governing  equa¬ 
tion  in  an  exponential  form,  logarithmic  interpolation  was  found  to 
produce  better  accuracy  than  linear  interpolation.  Sample  look-up 
data  that  were  generated  using  the  sub-grid  scale  model  is  shown 
in  Table  2. 

The  governing  conservation  equations  of  mass,  momentum  and 
current,  described  in  Section  2.2,  were  solved  using  a  conserva¬ 
tive  finite-volume  technique  [59].  The  SIMPLE  algorithm  [59]  was 
used  to  address  pressure-velocity  coupling  in  the  Navier-Stokes 
equations.  The  2D  model,  shown  in  Fig.  1,  was  used  for  simu¬ 
lations.  All  simulations  were  performed  on  a  uniform  grid  with 
50  cells  in  the  axial  direction  and  150  cells  in  the  cross-flow  (y) 
direction.  Nominally,  40  cells  were  used  across  each  flow  chan¬ 
nel,  10  across  each  GDL,  5  across  each  active  catalyst  layer,  and  40 
across  the  membrane.  This  particular  mesh  size  was  chosen  after  a 
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Simulation  parameters  and  values  of  key  properties  at  323  K  and  1  /I  atm. 
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Model  parameter 

Value/method  of  calculation 

Gas  channel  length 

7.112cm 

Gas  channel  width 

0.762  mm 

Diffusion  layer  width 

0.254  mm 

Membrane  width 

0.175  mm 

Membrane  permeability 

1.8  x  10-18  m2 

Diffuser  and  catalyst  layer  permeability 

1.76  x  10-11  m2 

Membrane  porosity 

0.28 

Anode  and  cathode  diffuser  layer 

0.5 

porosity 

Tortuosity  for  Bruggemann  correlation 

1.5 

Air  side/fuel  side  Pressures 

1/1  atm. 

Relative  humidity  of  inlet  streams 

100% 

Air  side  inlet  N2/O2  molar  ratio 

79/21 

H2  stoichiometric  flow 

2.8  A  cm-2  equivalent 

O2  stoichiometric  flow 

3.0  A  cm-2  equivalent 

Transfer  coefficients  at  anode  (Tafel 

0.5 

constants) 

Concentration  dependence  at  anode 

0.5  (H2) 

Reference  current  density  at  anode 

1.6  x  1011  (Am-3)  (m3/kmol  H2)1/2 

Membrane  electrical  conductivity 

Springer  et  al.  [57] 

Diffuser  layer  electrical  conductivity 

100  irr1) 

rigorous  grid-independence  study,  the  details  of  which  can  be 
found  in  Kamarajugadda  and  Mazumder  [42]. 

Moist  hydrogen  and  moist  air  was  introduced  into  the  anode 
and  cathode  inlet,  respectively.  A  plug  velocity  profile  was  imposed, 
and  the  magnitude  of  the  velocity  was  calculated  from  a  prescribed 
equivalent  current  density.  For  all  simulations,  the  temperature 
was  assumed  to  be  50  °C  everywhere,  and  the  energy  equation 
was  not  solved.  Other  relevant  parameters  used  for  simulations  are 
reported  in  Table  3.  The  kinetic  constant,  j0,  for  the  electrochemical 
reactions  at  the  anode  was  determined  by  calibrating  the  current 
density  at  low  bias  voltage  against  experimental  data  reported  by 
Ticianelli  et  al.  [60,61]— a  practice  that  has  been  used  in  the  past 
[15,42].  All  simulations  were  run  with  a  prescribed  bias  voltage 
boundary  condition  rather  than  a  prescribed  current  boundary  con¬ 
dition. 

All  transport  properties  of  the  fluid,  namely  viscosity  and  binary 
diffusion  coefficients,  were  computed  using  the  Chapman-Enskog 
equations  of  kinetic  theory  [46,62].  The  Lennard-Jones  potentials, 
which  are  needed  as  inputs,  were  obtained  from  the  CHEMKIN 
database.  The  density  of  the  fluid  was  calculated  using  the  ideal  gas 
law.  The  solutions  were  deemed  to  be  converged  when  the  residu¬ 
als  of  each  of  the  equations  decreased  by  seven  orders  of  magnitude. 
In  order  to  compute  the  actual  cell  voltage  from  the  prescribed  bias 
voltage  (or  potential  loss),  it  is  necessary  to  know  the  open  circuit 
potential.  As  done  by  previous  researchers  [15,63],  rather  than  use 
the  Nernst  potential,  an  empirical  correlation  [64]  was  used  for  the 
open  circuit  potential: 

0OC  =  0.00251  +  0.2329  (35) 

This  correlation  results  in  an  open  circuit  voltage  of  1.04  V  at 
50  °C.  Baseline  values  for  cathode  model  parameters  and  properties 
are  listed  in  Table  1. 

Fig.  7  shows  the  predicted  performance  with  cathodes  com¬ 
prised  of  agglomerates  of  intersecting  spheres  of  equal  radii.  The 
overlap  parameter  chosen  in  all  of  these  studies  is  §  =  0,  which  rep¬ 
resents  the  case  where  the  center  of  the  smaller  sphere  is  on  the 
surface  of  the  larger  sphere.  All  other  cathode  parameters  are  the 
same.  For  comparison,  the  predicted  performance  for  agglomerates 
comprised  of  isolated  spheres  is  also  shown.  The  results  indicate 
that  when  the  agglomerate  is  small  in  size,  the  exact  shape  does 
not  significantly  affect  the  overall  performance  of  the  PEMFC.  Fur¬ 
thermore,  the  shape  of  the  ionomer  (or  polymer)  coating  around  the 
nucleus  is  also  found  to  be  irrelevant,  as  also  corroborated  by  the 


Fig.  7.  Overall  PEMFC  performance  predicted  with  cathodes  comprised  of  agglom¬ 
erates  of  various  size  and  shape:  (a)  spheres  of  200  nm  radius,  (b)  spheres  of  600  nm 
radius,  and  (c)  spheres  of  1000  nm  radius. 


data  shown  in  Table  2.  In  this  case,  a  coating  thickness  of  <5/n  =  0.06 
was  considered.  When  the  agglomerate  is  large  [as  in  Fig.  7(c)],  both 
the  shape  of  the  nucleus  and  that  of  the  coating  plays  a  significant 
role  in  the  performance  at  high  current  density  (mass  transport  lim¬ 
ited  regime).  In  all  cases,  the  intersecting  spheres  are  found  to  result 
in  better  performance  than  a  single  equivalent  sphere,  in  keeping 
with  the  results  at  the  sub-grid  scale  (Section  3.2).  Therefore,  it  may 
be  concluded  that  the  better  performance  for  intersecting  spheres 
is  a  direct  result  of  increased  current  generation  at  the  agglom¬ 
erate  scale  due  to  increased  surface-area-to-volume  ratio  of  the 
agglomerate. 


338 


5.  Kamarajugadda,  S.  Mazumder  /  Journal  of  Power  Sources  208  (2012)  328-339 


Fig.  8.  Overall  PEMFC  performance  predicted  for  agglomerates  comprised  of 
1000-600  nm  radii  intersecting  spheres  with  £  =  1,  and  thin  and  capsule  coatings. 

Fig.  8  shows  the  predicted  performance  for  cathodes  comprised 
of  agglomerates  of  intersecting  sphere  of  unequal  radii— in  this  case 
1 000  nm  and  600  nm.  An  overlap  parameter  of  §  =  1  is  chosen,  which 
represents  the  nuclei  of  the  two  spheres  just  touching  each  other. 
At  the  agglomerate  scale,  this  particular  case  demonstrated  dra¬ 
matic  improvement  in  the  current  generated  when  compared  to  an 
equivalent  sphere,  as  has  been  depicted  in  Fig.  6.  Fig.  8  shows  that 
the  overall  performance  of  the  PEMFC  is  significantly  affected  by 
the  shape  of  the  agglomerate  as  well  as  the  type  of  coating  used.  As 
in  the  sub-grid  scale,  the  intersecting  spheres  result  in  better  per¬ 
formance  than  an  equivalent  sphere  of  the  same  volume,  primarily 
due  to  an  increase  in  the  surface-area-to-volume  ratio.  The  thin 
coating  produces  significantly  better  performance  than  a  capsule 
coating,  particularly  in  the  high  current  density  regime.  This  stems 
from  two  issues.  First,  a  thin  coating  poses  less  resistance  than  a 
capsule  coating  to  the  diffusion  of  dissolved  oxygen  to  the  reaction 
sites.  Secondly,  for  the  case  of  a  capsule  coating,  the  volume  fraction 
of  Nation  in  the  catalyst  layer  is  significantly  larger  than  in  the  case 
of  a  thin  coating.  Correspondingly,  the  volume  fraction  of  the  pores 
(average  porosity,  eCat)  is  significantly  lower  when  a  capsule  coat¬ 
ing  is  used.  For  example,  for  cathodes  comprised  of  agglomerates 
of  1 000-600  nm  intersecting  spheres,  the  porosity  was  calculated 
to  be  0.19  with  thin  coating  and  only  0.07  with  capsule  coating. 
This  poses  additional  mass  transport  resistance  to  the  diffusion  of 
oxygen  within  the  pores  of  the  cathode. 

4.  Summary  and  conclusions 

With  the  advent  of  next-generation  techniques  for  micro-/nano- 
fabrication,  it  will  be  possible  to  engineer  fuel  cell  catalyst  layers 
(electrodes)  to  desired  compositions  and  structures.  If  a  model  is 
available  that  is  general  enough  to  answer  the  question,  even  qual¬ 
itatively,  as  to  what  structure  and  composition  is  best  for  optimum 
performance,  its  impact  on  the  advancement  of  fuel  cell  technol¬ 
ogy  could  be  unprecedented.  The  model  developed  in  this  study  is 
a  logical  extension  of  the  popular  isolated-sphere  flooded  agglom¬ 
erate  concept,  and  takes  into  account  shape  effects  by  considering 
intersecting  spheres  of  unequal  radii.  Although  intersecting  spheres 
are  investigated  in  this  study  based  on  what  is  observed  in  micro¬ 
graphs  of  catalyst  layers,  the  model  is  applicable  to  agglomerates 
of  any  shape.  The  challenge  is  not  so  much  of  numerical  solution,  as 
of  parameterizing  the  shape  and  relating  it  to  an  actual  microstruc¬ 
ture.  Ultimately,  the  chosen  parameters  need  to  appropriately 
describe  the  real  microstructure,  or  some  statistical  representation 
thereof.  The  new  generalized  sub-grid  scale  flooded  agglomerate 


model  requires  numerical  solution  of  the  governing  reaction- 
diffusion  equation.  This  was  accomplished  using  the  unstructured 
finite-volume  method.  The  model  was  first  successfully  validated 
against  analytical  solutions  for  a  single  spherical  agglomerate.  Fol¬ 
lowing  the  validation  phase,  the  effects  of  agglomerate  shape  and 
size  on  the  overall  performance  (polarization  behavior)  of  a  poly¬ 
mer  electrolyte  membrane  fuel  cell  were  explored.  Based  on  the 
results  of  these  studies,  the  following  two  major  conclusions  may 
be  drawn: 

1.  Cathodes  comprised  of  agglomerates  of  intersecting  spheres 
result  in  a  better  overall  fuel  cell  performance  than  cathodes 
comprised  of  agglomerates  of  individual  spheres  of  the  same 
equivalent  volume.  This  is  a  manifestation  of  the  improvement 
in  the  oxygen  reduction  rate  (or  current  generation  rate)  at 
the  agglomerate  scale  due  to  increased  surface-area-to-volume 
ratio. 

2.  For  small  agglomerate  sizes,  the  spherical  flooded  agglomerate 
model  with  isolated  equivalent  spheres  (as  used  in  all  previ¬ 
ous  studies  using  the  flooded  agglomerate  model)  appears  to 
be  an  acceptable  approximation  even  when  the  agglomerates 
are  not  strictly  spherical  in  shape.  In  other  words,  shape  effects 
are  insignificant  when  the  agglomerates  are  small.  Flowever,  for 
larger  agglomerates  (>600  nm  radius),  shape  effects  are  signifi¬ 
cant,  and  a  model  that  accounts  for  the  shape  of  the  agglomerate, 
as  the  one  proposed  here,  is  necessary  for  accurate  prediction  of 
PEMFC  performance. 

In  the  future,  the  study  will  be  extended  to  include  the  effect 
of  liquid  water  and  non-isothermal  effects.  Flowever,  the  authors 
hope  that  this  study  has  been  able  to  provide  some  directions  for 
future  improvements  in  modeling  the  cathode  of  a  PEMFC  without 
actual  microstructure  reconstruction. 
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